Spatial Distributions of Observables in Systems under Thermal Gradients 
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Departures of observables from their thermal equilibrium expectation values are studied under 
heat flow in steady-state non-equilibrium environments. The relation between the spatial and tem- 
perature dependence of these non-equilibrium behaviors and the underlying statistical properties 
are clarified from general considerations. The predictions are then confirmed in direct numerical 
simulations within the FPU-/3 model. Non-equilibrium momentum distribution functions are also 
examined and characterized through their cumulants and the properties of higher order cumulants 
are discussed. 
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I. INTRODUCTION 

In studies of non-equilibrium physics, especially those of steady states, local equilibrium is most often invoked and 
this assumption simplifies calculations through the use of equilibrium statistical mechanics and thermodynamics^. 
The local equilibrium assumption allows the use of the equilibrium distribution function to compute observables. 
If local equilibrium conditions are not assumed, very little can be computed analytically and even the definition 
of temperature is no longer unique 0, 0- Efforts have been made to quantify the goodness of local equilibrium 
assumptions or how trans por t coefficients differ from their linear response values, though only few quantitative studies 
exist 1 S @, 0, l lift . Without the knowledge of the non-equilibrium steady-state distribution, theoretical 
development becomes quite restrictive. We explore how observables depart from their equilibrium expectation values 
within a given non-equilibrium steady-state, specifically focusing on the spatial dependence of the non-equilibrium 
expectation values within a given system and their local temperature dependence. To make this concrete, heat flow in 
the FPU (3— model is simulated to test the predictions. We further quantitatively examine the relationship between 
the momentum cumulants and the distribution and find that the lower order cumulants characterize the distribution 
quite well. 

For systems in thermal gradients, it is natural to consider how an observable O in the non-equilibrium steady state 
departs from its equilibrium value, denoted O eq . The normalized deviation from equilibrium, when O eq ^ 0, can be 
expanded as 
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When O eq = 0, as is the case for higher order momentum cumulants, one can normalize by an observable which has 
the same dimensions. When local equilibrium is no longer valid, in general, no unique definition of temperature exists 
and a choice needs to be made. This definition of non-equilibrium temperature can be thought of as a choice of a 
coordinate system, on which the physics behavior of the system will not depend. If we assume analyticity in VT, 
the deviations So can be expanded in even powers as above. We shall see below that this expansion is adequate for 
describing the properties of the system. 

The heat flow, J, is the flow of energy and can be unambiguously defined in Hamiltonian systems. Near equilibrium, 
it satisfies Fourier's law locally as J = — kVT(i), where k is the thermal conductivity, T(x) is the temperature profile 
inside, and x is the position inside the system. Fourier's law can be used in Q to re-express the local departures 
from equilibrium in terms of the temperature profile T{x), or equivalently the position x once the coefficients C, C 
are known since J does not depend on x; 
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We note that Fourier's law itself receives non-equilibrium corrections which is why the coefficient of 0(J 4 ) term 
in the expansion J5J differs from that of QJ. I n the following, the objectives will be to make the formula more explicit 
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and understand its physical properties under rather general assumptions. This relation, together with k(T) (and 
consequently T(x)) provides the basis for defining how non-equilibrium observables vary inside a finite system both 
near and far from global thermal equilibrium. 



II. THE FPU MODEL AND TEMPERATURE PROFILES 



The results we present here are derived from general considerations and we develop them in conjunction with a 
model in which they can be explicitly analyzed. We study the FPU (3 Hamiltonian, defined generally in the form 
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We use the FPU model since its physical properties are of wide interest f [T2I Il3l [Til Il5j and references therein). Also 
as the model is well studied, we can understand the physical properties we find within a larger physics context. Under 
the rescaling pk = p' k L0 2 Vm 3 , qk — q k w*/m, we obtain the conventional form of the FPU (5 model, 
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where Hp = H / (m 2 ui 1 ). We note that in finite temperature simulations, changing the temperature is equivalent to 
changing the coupling (3. Under the additional rescaling p' k — pk/ \f]3, q' k = qk/ VP, one obtains a unique, dimensionless, 
Hamiltonian H = Hp—i — (3Hp, which we shall use without any loss of generality. Since p 2 — f3p' 2 , the temperatures 
in the two formulations H and Hp are related by T = f3T'. 

In this work, we study the non-equilibrium steady state physics of the theory under thermal gradients, making 
use of non-equilibrium states constructed numerically. (For general discussion, see, for instance, 03E3.) The 
model is thermostatted at the boundaries k = 0, L at various temperatures T-j 3 , T§ , using the generalized versions 
of Nose-Hoover thermostats as detailed in . These additional thermostat degrees of freedom are added only at 
the boundaries and the degrees of freedom inside the system (0 < k < L) are exclusively those of the Hamiltonian 
(@J. By numerically integrating the equations of motion of the whole system (including those of the thermostats), 
we obtain the behavior of physical observables in the non-equilibrium steady state by averaging over time, in the 
standard manner [l7j. The local temperature at site k is defined as T k = (p 2 )- In this work, we study the physics 
inside the system, away from the boundaries by much more than the mean free path of the system^]. The sensitivity 
of the results to the manner in which we apply the boundary conditions — including both the number of thermostats 
and the strength of the couplings — have been examined to ensure that physics results below remain independent of 
their implementation. (The only exceptions are the boundary jumps in temperature which we discuss below.) The 
numerical integrations were performed using the fourth order Runge-Kutta routines with time steps of 0.005 ~ 0.02 
for 10 7 ~ 10 10 time steps. The equilibrium properties have been readily verified with this methodjlj, I3- 

In Fig. ^ some examples of temperature profiles for the FPU theory are shown. Generically, there are temperature 
jumps just inside the boundaries with smooth temperature variations within. The boundary jumps become larger 
as one moves away from global equilibrium. The jumps are dynamical in the sense that they depend on the model, 
the transport coefficient, heat flow, as well as the type of boundary conditions employed. The temperatures at the 
boundaries are at the thermostat temperatures to high degree of precision. For instance, in the examples of Fig. ^ 
the boundary temperatures are equal to the prescribed thermostat temperatures to within few in 10 5 relatively. 

From temperature profiles and heat flow calculations, Fourier's law can be verified to hold up to corrections of the 
form l[TJl. and the thermal conductivity, K, can be obtained for a given temperature and system size. In the 1-d FPU 
model, K depends on the system size L and does not display bulk behavior 13J. k is also dependent on the temperature 
in a known manner [l4|. Generally, in cases where we have a one dimensional temperature gradient, the temperature 
profiles can be obtained by integrating Fourier's law as long as we are not too far from equilibrium [111 Il8l Il9| : 



K (T) dT = -Jx, J = -p k [(%+i - qk) + (qk+i - qkY 



(5) 



x is the continuum extrapolation of the discrete lattice index k. We note here that J is a constant within the system 
for a given set of temperature boundary conditions since there are no heat sinks or sources inside. T\ in the integral 
is the temperature extrapolated to the boundary and is explained below. 
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FIG. 1: Some examples of temperature profiles for the FPU model with L — 128. The thermostat temperatures at the 
boundaries are (T?,T$) = (0.88, 16.72), (2.4, 15.2), (4.4, 13.2), (6.6, 11.0) for the four thermal profiles. The profiles predicted 
from Eq. J7J are indicated by x and agree well with the results from the numerical simulations. 



In many situations, the temperature dependence of the thermal conductivity, within some temperature range, can 
be well described by 



k(T) = cT- 



(0) 



While this power law may not hold globally in T, it is often the case that it is sufficient for the region of interest, 
which is the case here. In such a situation, the temperature profile can be explicitly computed from (|SJ to be|l8|| 
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Here, T± t 2 denote the boundary temperatures obtained by extrapolating the temperature profile inside the system 
and differs from the thermostat temperatures T® 2 by the boundary temperature jumps. From (jSJ) and JSJl, the 
temperatures T\ t 2 are found to obey a relation 
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To understand the temperature profile of the whole system, we further need an understanding of the temperature 
jumps at the boundaries [2(J. Similar boundary slips have been seen in sheared systems and these effects have been 
known for a long time in real systems. To leading order, the temperature jumps can be described by (with n being 
the normal to the boundary) 
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Here A is the mean free path of the excitations, which for the FPU lattice model, is essentially the k(T) (up to a 
constant factor of order one) due to kinetic theory arguments a reflects the efficacy of the boundary conditions. 
The last relation is obtained by using Fourier's Law and (JSJ). The jumps on the hot and cold side are the same provided 
the system is reasonably close to equilibrium. The jumps at the boundaries and the temperature profile within lj?Jl 
describe the temperature profile of the complete system. The predicted values for the temperature profiles are plotted 
in Fig. ^ at a number points inside the systems (x symbols) away from the boundaries and are seen to be consistent 
with the simulation results. The thermal conductivity is roughly constant with respect to the temperature in this 
region so that 7 = was used in the profile calculations. This demonstrates that all aspects of the non-equilibrium 
temperature profile can be quantitatively captured through Q and @, irrespective of whether n(T) is a power law in 
temperature for all T or not. With this understanding of T(x) we can now turn to the question of general observables. 
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TABLE I: Non-Equilibrium coefficients C 4 = {fJ.T s )L a for ((p 4 
model and the 4 theory in d — 1 ~ 3-dimensions. The value c 
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/T 2 (c/. Eq. (HJjjni)- The results are shown for the FPU 
s is extracted from fitting to several temperatures. 



III. SPATIAL DEPENDENCE OF CUMULANTS IN THE NON-EQUILIBRIUM STEADY STATE 

In non-equilibrium steady states, physical observables show deviations from their equilibrium values reflecting the 
lack of local equilibrium in the system. The behavior of the observables have been seen to be well described by 
on average, at least in some cases |ll|. Here, we now would like to investigate a more detailed issue — whether these 
properties can be used to understand the nature of the spatial profiles of these observables in a given non-equilibrium 
situation. We will assume that within some range of T and L that we can represent the expansion coefficients in 
as 

C = PoT s °L a ° (10) 

The behavior of Co with respect to T, L clearly must depend on the dynamics of the theory and is not expected to 
be generic. 

To study the spatial distribution of physical observables in non-equilibrium, we make use of J2J which describes 
how the observables should behave in non-equilibrium locally in space, given the thermal conductivity. Using this 
property and ©, we obtain to leading order that observables will deviate from their local equilibrium values as 

So = Co ( JT W y ~ 1 ) 2 = aoT(x)^-» + > (11) 



Here a® is defined through this equation and should be proportional to J 2 . This implicitly contains the spatial 
distribution since the temperature profile is known and can be understood as in (JJJ. 

While these arguments apply to any physical observable in the system, we choose to study cumulants of momenta, p, 
mainly for the following reasons; conceptual and practical. There seems to be no universal rigorous definition of local 
equilibrium, yet the concept in the least seems to include a unique meaning for temperature, which in this case would 
lead to the Maxwellian distribution for p. To put another way, when the momentum distribution is not Maxwellian, 
we can choose different definitions of the temperature based on the various moments of p|2ild|- The cumulants of the 
momentum distribution provide insight into how the physical properties of a non-equilibrium system deviates from 
those of local equilibrium. The cumulants are well defined local variables and their values in local equilibrium are 
known precisely. The low order cumulants are defined as 

«p 2 » = (p 2 ), «p 4 » = (p 4 > - Hp 2 ) 2 , «p 6 » = <p 6 } - iHp 2 )(p 4 ) + 30( P 2 ) 3 , . . . (12) 

where, in equilibrium, 

((p 2 )) eq = T, «p"» e9 = (n + 2) (13) 

This property is also of practical importance. Since the deviations we compute can be small, it is desirable to use 
observables whose local equilibrium values are known exactly. In this case in thermal equilibrium, O eq = 0, so we 
use So = ((p 2n ))/T n . We list the coefficient for the case O = ((p 4 )) in Tabledfor the FPU /3-model as well as (j) 4 
theory [l8j for comparison. 

Let us investigate how well (|11|) describes the spatial distribution of ((p 4 ))/T 2 . We find 
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FIG. 2: (left) Spatial dependence of the rescaled 4-th momentum cumulant, ({p 4 )) /T 2 for the four systems in Fig. Q Larger 
cumulant values are seen for larger boundary temperature differences, (right) Temperature dependence of {{p 4 )) /T 2 for the 
same systems. In both panels, the predictions 111 U are indicated by the dashes. 
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FIG. 3: J dependence of the non-equilibrium expansion coefficient, an, for various boundary conditions, (Ti,T§) and system 
sizes L. The dashed line is 3.72 J 2 and the ~ J 2 behavior of the coefficient can be clearly seen, as predicted from theory. Each 
data point represents a particular temperature boundary condition for L — 32 (x) , L — 64 (□) and L = 128 (Q) systems. 



S4 is the temperature dependence of the coefficient C4 which is reflected in Table [I] To understand the validity of 
the prediction Eq. Q14[). fits were made with just one parameter 0,4 for the whole profile. We find that this describes 
the situation quite well, as seen in the examples of Fig. |2J where the predictions are denoted by dashes. In these 
figures, we have compared the fits with the spatial as well temperature dependence of ((p 4 )) for the four systems 
shown in Fig. ^ In this temperature range, temperature dependence of the thermal conductivity is weak so we used 
7 = and S4 = —0.14 extracted from the data in Table [I] Similar results were found for different temperature 
boundary conditions and for different L. To further verify the underlying physics, we study the J dependence of the 
coefficient 04. The behavior for various systems, including the four systems in Fig. ^ are shown in Fig. [3] Each data 
point represents a system with a particular size and temperature boundary conditions. The central temperature is 
around T = 8.8 and is kept fixed. The observed behavior is clearly well described by 04 ~ const. x J 2 . The coefficient 
04 seems //-independent and this can roughly be understood since c 2 grows in L in a manner similar to C4. We 
have in addition systematically studied the results to see if we can discern the contribution of higher order terms in 
the expansions 0,1© (of order J 4 and higher) but have found no consistent evidence for them. In other physical 
situations, non-analytic behavior seems to have been seen in some cases pll l22| . 

While the logic seems to work for the lowest non-trivial order cumulant, ((p 4 )), we find it instructive to analyze if 
it works at higher orders. In this direction, we have analyzed the next non-trivial order ((p 6 )) and have found that 
its behavior is quite consistent with physics of (|llf> . as was the case of ((p 4 )), in all the systems we have studied. In 
Dractice. higher order cumulants are more Drone to errors and the conrDutations are more difficult. The results for 



6 




100 120 




FIG. 4: Spatial dependence (left) and temperature dependence (right) of ({f> 6 ))/T 3 for the four systems in Fig. Q Larger 
cumulant values are seen for larger boundary temperature differences. Predictions are shown with dashes. 
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FIG. 5: J dependence of the coefficient ag for various boundary conditions, (Ti 1 ,!^) an d system sizes L. The dashed line 
denotes 156 J 2 . ~ J 2 dependence of ae is evident, in agreement with the predictions. Each data point represents a particular 
temperature boundary condition for system sizes L = 32 (x), L = 64 (□) and L = 128 (O)i as m Fig- El 



the same four systems in Fig.^are shown in Fig.|I] As in the ((p 4 )} case, the coefficient a?, shows J 2 behavior within 
error, as it should, ag shows a weak L dependence, as we would generically expect. A common value of sq = —1.6 
was adopted for all the data in Fig. 0] and Fig. |SJ What is evident is that the spatial behavior of non-equilibrium 
observables can be explicitly related to transport and other physical properties of the system using rather general 
considerations. From the cumulants we now consider what can be said about the full momentum distribution function. 



IV. CUMULANTS AND THE DISTRIBUTION 



The cumulants are quantitative indicators of the non-Maxwellian nature of the momentum distribution or the 
violations of local equilibrium. All the cumulants are non-zero unless the system is in local equilibrium, in which 
case only the linear and quadratic cumulants are non-zero. There are very few problems where cumulants can all be 
computed analytically and it becomes numerically intractable to compute them as we go to higher orders. It is then of 
interest to see how well the lower order cumulants characterize the distribution. The cumulants are properties of the 
distribution function, which has an infinite number of degrees of freedom. A priori, there is no reason to assume that 
the lower order cumulants characterize the distribution. In order to clarify this issue, first note that the distribution 
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function f(p) and the cumulants are related explicitly through the generating function as 



Here, in the last equality, the symmetry under p «-► — p was used, which leads to {(p 2n+1 )) — 0. We see from 
this equation that given all the cumulants (or equivalently, moments), we may recover the distribution function by 
performing an inverse Fourier transform. However, in practice, not all the cumulants are available. 

Intuitively, we expect the lower order cumulants to be the leading order results with higher order cumulants becoming 
more important as we move further away from equilibrium. In Fig. |B| we plot the relative difference of the measured 
distribution f(p) to the thermal distribution, /o(p), for the distribution directly measured in the simulations and the 
distribution computed from the low order cumulants, ((p 2 ' 4 ' 6 )). The comparisons are performed for the four systems 
in Fig. at a point in the middle of the system. From these graphs, we observe the following: (a) The agreement 
between the distribution computed from lower order cumulants and the distribution is quite good in all cases; (b) 
the relative deviation from the thermal distribution is larger as we move away from equilibrium (larger AT/T), as 
expected; (c) the small discrepancy between the computed distribution and the measured one seems to be larger 
for larger AT/T; (d) the deviation from the thermal distribution becomes more noisy for smaller AT/T, since the 
deviation itself is smaller and the relative error is larger. We mention here that strictly speaking, the distributions can 
have different behavior, such as long tails, beyond the region we have investigated. However, these tails would have 
to be quite small since the distributions decay as exp(— p 2 /(2T)) and the agreement is good up to reasonably large 
p, as seen in Fig. [HI We have examined numerous systems for different T and L and found similar good agreement. 
Therefore, we see that the lower order cumulants provide good physical observables that quantitatively describe the 
deviations of the systems from local equilibrium, at least in the FPU model. 

It is possible to examine the characteristics of the higher order cumulants. It should be noted that unlike the even 
moments (p 2n ), even cumulants, ((p 2 ™)), need not be positive and in general will not be. So to study the general 
trend of the cumulants for higher order, we examine the magnitude of the cumulants. In Fig. (left), we show the 
behavior of the cumulants up to 20-th order for the same four systems in Fig. specifically for the point at which 
the momentum distributions in Fig. were computed. Only data points with reasonable error are shown and an 
explanation of the relevant errors is given below. We see an increase in the magnitude with the order is roughly 
exponential. This growth is far milder than the (2n)\ seen in (|15|l . 

The behavior of the higher order cumulants is of some import and we briefly explain semi-quantitatively why they 
are difficult to obtain. The difficulty lies mainly in the statistical error in the simulations. This can be estimated from 
the number of samples for computing the expectation values as 

^ ~ 4= (16) 

where A denotes the error and N is the total number of samples or the number of time steps in the simulation. Note 
that ((p™)) = (p n ) + ... so that an error estimate for the moment should suffice as the error estimate for the cumulant. 
An adequate value for the moment can be obtained in equilibrium, 

^ -• (n-l)U (17) 



T n/2 



Combining these relations, we find the statistical error for the cumulants which increases rapidly for higher order 
cumulants. 

(n-l)Hn 



These estimates for the error also apply to the equilibrium situation. In contrast to the non-equilibrium cumulants, 
the equilibrium cumulants should vanish, with the exception of ((p 2 )). As the measured values will converge to zero, 
at any given time-step in the simulation, their values will be generically non-zero. In Fig. (right), we compare 
the equilibrium cumulants, in the middle of the system to the above error estimates. It can be seen that the rough 
estimate (|18ll seems to be consistent with the results. As one samples more [N increases), these will tend to zero. 
However, for a finite sample size, this is found to explain the order of the uncertainty. 

With N = 10 9 time-steps — which we used for the values in Fig. — for 8 and 10-th order cumulants, the errors 
are 0.03 and 0.3. As we can see from Fig. (left), this means that we can obtain up to the 8 or 10-th cumulant with 
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FIG. 6: The relative deviation of the distribution from the Maxwell distribution for the four systems in Fig. Distribution 
obtained from the cumulants {(p 4 )}, ((p 6 )) (dashed) are compared with the measured distributions (solid). The agreement is 
excellent. AT/T denotes the boundary temperature difference over the average temperature and is an indication of how far 
the system is from equilibrium. 



to equilibrium. These error estimates are quite consistent with the estimates we obtain from the statistical properties 
of the simulations. These errors can be overcome with higher statistics which quickly becomes unrealistic for higher 
order. We have analyzed systems with various other temperature boundary conditions and L and have found the 
increasing behavior of the cumulants seen in Fig. □ (left) to be quite generic. 



V. SUMMARY AND DISCUSSIONS 



The spatial distribution of cumulants in non-equilibrium steady states under thermal gradients were predicted from 
general considerations and tested in the the FPU model. The understanding of the temperature profile for a given 
non-equilibrium steady state, combined with the deviations of physical observables from their equilibrium values, can 
be used to develop a consistent description of the spatial distribution of observables. In principle, the behavior of 
observables probably have higher order corrections in the non-equilibrium nature of the system, which in this case is 
VT, but higher order effects could not be separated within the current numerical simulation results. 

We quantitatively analyzed the relation between the momentum cumulants and the distribution in the non- 
equilibrium steady state. It was found that the lower order cumulants characterize the difference of the non-equilibrium 
distribution from the one in local equilibrium quite well. Understanding and characterizing the properties of the dis- 
tribution is of manifest importance since the distribution function for physical variables allows us to compute any 
observable constructed from these variables. To understand the properties of any local variable in the non-equilibrium 
state, the physical properties of the coordinate variables also need to to be clarified. 

A comment is perhaps in order: lack of local equilibrium behavior can in some cases be attributed to the lack of 
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FIG. 7: (left) Higher order cumulants, \{{p"))\/T n/2 (n < 20) for the four systems in Fig. □ (T?,T 2 °) = (0.88,16.72), (x), 
(2.4, 15.2), (□), (4.4, 13.2), (O) an d (6.6, 11.0), (A). Only points with reasonably small error are shown. The dashed line is 
5.0 x 10 -6 exp(1.5n) drawn for comparison, (right) The equilibrium cumulants for L — 16 (x), L = 32 (□), L = 64 (A) and 
L = 128 (O) compared to the rough estimate, Eq. I)18|l (dashes). The cumulants ((p 2 )) were measured at the middle of the 
system with number of samples iV = 10 9 atT = 8.8. 



coarse graining 4]. Heuristically speaking, if one does not have a large number of degrees of freedom, one cannot see 
the equilibrium behavior. This is a different phenomenon from the case at hand, since the effective number of degrees 
of freedom is the number of samples in the ensemble average which is taken in the time averaging procedure. This 
number is huge. In fact, as is well known, in these types of ensembles, it makes perfect sense to talk even about the 
statistical mechanics of one spin degrees of freedom. This is also quite clear from our results; the deviations from local 
equilibrium seen in Fig. EI Fig- El Eq. J3J are of definite sign and no amount of averaging over space will make it zero. 
So coarse graining will not average out the violations of local equilibrium seen above. Also, the non-local equilibrium 
properties found in this paper pertain to systems in the non-equilibrium steady state and therefore are not transient. 

We have also performed similar analyses of spatial distributions on the (f) A model. The physical properties of the 
model are different from those of FPU model and we intend to report on this in the near future. 
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